Algoritmos de Aprendizaje Automático con Datos Espaciales
Ing Agr. (Dr.) Mariano Córdoba
Ing Agr. Franca Giannini Kurina
CONICET, UNC-Argentina
9 y 10 de octubre, 2019
Módulo 2:
Algoritmos de Aprendizaje Automático.
Random Forest.
Ajuste y validación de modelos predictivos.
Incorporación de la información espacial: Random Forest Kriging.
Ejemplos de aplicación.
Estadística Clásica vs. Geoestadística
Considera los datos como realizaciones independientes de una misma variable aleatoria (VA) (siguen la misma distribución)
Se busca estimar los parámetros de esta distribución
Cuando los datos están ubicados en el espacio geográfico estas hipótesis son raramente aceptables
Observaciones próximas en el espacio tienen valores “parecidos”, mientras que aquellas que están más alejadas entre sí tienen una menor relación
Variable Regionalizada
Un conjunto de valores de una variable aleatoria constituyen una función aleatoria o proceso estocástico.
La realización de una función aleatoria se conoce como variables regionalizadas
El fenómeno regionalizado se extienden en el espacio y presenta una cierta continuidad
Las variables regionalizadas miden propiedades o atributos relacionados con la variación del mismo en el espacio
Algunos ejemplos
Porosidad de la roca en un reservorio de petróleo o en un acuífero
Concentración de un elemento contaminante en la atmósfera o en el suelo
La elevación en un punto del espacio geográfico
El pH y la concentración en nutrientes medidas sobre una muestra de suelo
El número de árboles y su diámetro promedio en áreas de observación de un bosque
Valor unitario de la tierra
Tendencia respecto a coordenadas
Distribución de VUT
## [1] 1.795173
Visualización en el plano
Visualización en el plano
Proceso Aleatorio
Estos procesos están físicamente determinados, pero sus interacciones son tan complejas que la variación parece ser aleatoria
Dada la complejidad de la interacción, una solución determinística puede distorsionar la inferencia estadística
La solución propuesta por Matheron (1965) fue tratar la variación como si fuera aleatoria
Estacionariedad
La noción de estacionariedad sostiene la geoestadística y nos permite asumir que existe el mismo grado de variación de un lugar a otro.
El proceso aleatorio puede ser representado por el modelo \[Z(x)=µ + ε(x)\]
donde \(\mu\) es la media del proceso y \(ε(x)\) es una cantidad aleatoria con una media de cero y una covarianza, \(C(h)\), dada por
\[C(h)= E[ε(x) + ε (x+h)]\]
\[C(h)= E[(Z(x) -µ) (Z(x+h) - \mu) ]= E[Z(x) Z(x+h) - \mu^2]\]
Estacionariedad
Si la media no es contaste, entonces la covarianza no puede existir
En ese caso invocamos la hipótesis de estacionariedad débil o intrinseca propuesta por Materon en la que las diferencias esperadas son cero,
\[E[Z(x)-Z(x+h)]=0\]
La co-varianza se sustituye por la mitad de las diferencias, la semivarianza:
\[\gamma(h)={\frac{1}{2}}var[Z(x)-Z(x+h)]={\frac{1}{2}}E[(Z(x)-Z(x+h))^2]\]
Al igual que la covarianza, la semivarianza depende de \(h\) y sólo de \(h\), y como la función de \(h\) es el variograma \(\gamma(h)\)
Estacionariedad
Un proceso estacionario a una escala podría, a otra escala, mostrar tendencia, es decir, un componente sistemático
el modelo para esa escala sería:
\[Z(x)=\mu(x)+e(x)\]
En la que \(\mu(x)\) es un término de tendencia determinista que sustituye a la media constante \(\mu\)
El variograma es calculado a partir de los residuos
\[\gamma(h)={\frac{1}{2}}E[(e(x)-e(x+h))^2]\]
Semivarianza
Función básica que describe la variabilidad espacial de un fenómeno
Responde a la siguiente pregunta:
¿Cómo de parecidos son los puntos en el espacio a medida que estos se encuentran más alejados?
Mientras más similares sean las observaciones, menor la semivarianza
Semivarianza
\[\gamma^*(h)={\frac{1}{2N(h)}}\sum_{i=1}^{n}\big[(E(x)-E(x+h))^2\big]\]
\(\gamma(h)\) es el valor experimental de la semivarianza de datos separados por una distancia \(h\) \(N\) es el número de pares de puntos de muestreo separados por la distancia \(h\)
\(Z(x_i)\) es el valor de la variable \(Z\) en el punto muestreado \(x_i\)
\(Z(x_i + h)\) es el valor de la propiedad \(Z\) en el punto muestreado a una distancia \(x_i + h\)
Semivariograma experimental
El semivariograma es generalmente una función incremental. A distancias crecientes entre 2 puntos la varianza de su diferencia tiende a aumentar.
Se definen “ventanas” en torno a una unidad de muestreo a los fines de encontrar los puntos pares correspondientes a distancias y direcciones determinadas.
Semivariograma Teórico
Algunos de los modelos teóricos más usados para modelar semivariogramas asumen tendencias lineales cerca del origen (modelos exponenciales y esféricos) mientras que otros involucran un efecto cuadrático cerca del origen (modelo gaussiano).
Un incremento lento del semivariograma cerca del origen sugiere la presencia de un patrón de estructura espacial consistente.
Para un mismo rango y umbral, el semivariograma exponencial crece más rápidamente que el gaussiano, sugiriendo menor estructura espacial.
Ajuste de Semivariogramas
Modelo exponencial
## model psill range
## 1 Nug 35520.78 0.000
## 2 Exp 3514832.04 1022.817
Ajuste de Semivariogramas
Modleo esférico
## model psill range
## 1 Nug 212383 0.000
## 2 Sph 2973759 2039.737
Visualización conjunta
vgLine <- rbind(cbind(variogramLine(mod_exp, maxdist = max(semi_exp$dist)), id ="Exponencial"),
cbind(variogramLine(mod_esf, maxdist = max(semi_exp$dist)), id ="Esférico")
)
grafico <- ggplot(semi_exp, aes(x = dist, y = gamma, color = id)) +
geom_line(data = vgLine) +
geom_point() +
labs(title ="Semivariograma experimental y teórico ajustado") +
xlab("Distancia") +
ylab("Semivarianza") +
scale_color_discrete(name = "Modelo", labels = c("Esférico", "Exponencial", "Experimental"))
grafico
Evaluación y ajuste automático
## [1] 268556207
## [1] 272019985
## [1] 268556210
Interpolación espacial
Poligonos de Thiessen (Voronoi)
Triangulación
Inverso de la distancia (IDW)
Regresión polinómica…
Kriging
Interpolación Kriging
El kriging es una técnica de estimación que proporciona el mejor estimador lineal insesgado (BLUE, Best Linear Unbiased Estimator)
Proporciona un error de estimación conocido como varianza de kriging que depende del modelo de semivariograma obtenido y de las localizaciones de los datos originales
Esto permite hacer un análisis sobre la calidad de las estimaciones
El método kriging realiza una ponderación, i.e. atribuir un peso a cada valor observado
Los pesos son calculados de manera que minimice la varianza de estimación resultante, teniendo en cuenta la característica espacial de los datos
Al minimizar la varianza de estimación se garantiza el uso óptimo de la información disponible
Interpolación Kriging
Supongamos que los valores de la variable aleatoria \(Z\), han sido muestreados en los puntos \(x_1,x_2,...,x_n\)
para \(N\) datos \(z(x_i)\), \(i=1,2,3,...,N\) y se quiere predecir \(Z\) en cualquier punto \(x_0\), mediante: \[\hat Z(x_0)\sum_{i=1}^{N}\lambda_{i} z(x_i)\] Donde \(\lambda\) son los pesos asignados para la ponderación
Interpolación Kriging
Para asegurar que la estimación sea insesgada los pesos son dados de manera que:
\[\sum_{i=1}^{N}\lambda_{i}=1\] \[E[\hat Z(x_0)-z(x_0))]=0\]
Siendo además un estimador de varianza mínima:
\[var[\hat Z(x_0)]=E[(\hat Z(x_0)-z(x_0))^2]= \sum_{i=1}^{N}\lambda_{i}\gamma(x_i-x_0)-\sum_{i=1}^{N}\sum_{j=1}^{N}\lambda_{i}\lambda_{j}\gamma(x_i-x_0)\]
Grilla de predicción regular
Grilla de predicción regular
Grilla de predicción por ejes de manzana
## [1] "g_gas" "g_agua" "g_cloaca" "g_perc_baldio"
## [5] "g_comercio" "g_sellos" "g_perc_edif" "d_bpriv"
## [9] "d_centro" "d_educ" "d_espverd" "d_ffcc"
## [13] "d_supcom" "d_indust" "d_ruta" "d_viasprin"
## [17] "d_viassec" "d_depre" "d_front" "x"
## [21] "y" "g_vial" "d_sube" "cod_baja"
## [25] "d_baja"
## [1] "VUT" "d_bpriv" "d_centro" "d_educ" "d_espverd"
## [6] "d_ffcc" "d_supcom" "d_indust" "d_ruta" "d_viasprin"
## [11] "d_viassec" "d_depre" "d_front" "g_comercio" "g_por_bald"
## [16] "g_porc_edi"
Grilla de predicción por ejes de manzana
Kriging
Simple -> variable estacionaria de media conocida
Ordinario -> variable estacionaria de media desconocida
Universal
Co-Kriging
Regression Kriging
Kriging con deriva externa
Kriging
Si el número de variables auxiliares es bajo y no están disponibles en todos los nodos de la grilla, se debe usar co-kriging
Si se dispone de información auxiliar en todos los nodos de la grilla y se correlaciona con la variable objetivo, se debe usar kriging con un modelo de tendencias, también llamado "deriva externa“.
- Kriging universal: se produce un cambio de la media cuando la tendencia se modela en función de coordenadas.
- Kriging con deriva externa: la deriva se define externamente a través de alguna variable auxiliar
Kriging en la grilla
## [using ordinary kriging]
Kriging en la grilla
Kriging en los ejes de manzana
## [using ordinary kriging]
Kriging en los ejes de manzana
Vemos todo junto!!
# EXtracción de predichos Y varianza de predicción
pred <- kriging_grilla$var1.pred
varpred <-kriging_grilla$var1.var
cord_X <- kriging_grilla$x
cord_Y <- kriging_grilla$y
tablapred <- cbind(cord_X,cord_Y,pred, varpred)
# Visualización interactiva
cols <-colorRampPalette(rev(brewer.pal(11, "Spectral"))) # Genero una paleta
muestra <- mapview(datos,zcol = "VUT", ceX ="VUT", layer.name ="Muestra", col.regions = cols, alpha=0)
pred_grilla <- rasterFromXYZ(tablapred[,1:3])
crs(pred_grilla) <-CRS("+init=epsg:22174")
mapapred_grilla <- mapview(pred_grilla,legend=T, col.regions = cols, layer.name ="Predichos Grilla",na.color = "transparent")
var_grilla <- rasterFromXYZ(tablapred[,c(1,2,4)])
crs(var_grilla) <-CRS("+init=epsg:22174")
mapavar_grilla <- mapview(var_grilla,legend=T, col.regions = cols, layer.name ="Varianza Predicci?n Grilla",na.color = "transparent")
crs(kriging_ejes) <-CRS("+init=epsg:22174")
mapapred_ejes <- mapview(kriging_ejes, zcol = "var1.pred", ceX ="var1.pred",col.regions = cols, layer.name ="Predicho Ejes", alpha=0)
mapavar_ejes <- mapview(kriging_ejes, zcol = "var1.var", ceX ="var1.var",col.regions = cols, layer.name ="Varianza Predicci?n Ejes", alpha=0)
muestra + mapapred_grilla + mapavar_grilla + mapapred_ejes + mapavar_ejes
Otras opciones de visualización
Otras opciones de visualización
Muchas gracias
Nos vemos mañana :)
